The invention relates to a computer for evaluating signals from nuclear magnetic resonance tomography and a nuclear resonance tomograph equipped with such a computer.
In medical research, it is desirable to obtain information with respect to the brain activity or, in a wider sense, information concerning the blood flow or the deoxyhemoglobin concentration changes in animal and human organs. The neural activation results in an increase of the blood flow in activated brain areas, which causes a decrease in the blood deoxyhemoglobin concentration. Deoxyhemoglobin (DOH) is a paramagnetic compound which accelerates the nuclear spin-dephasing. If the DOH concentration is lowered because of a brain activity which causes a blood flow, the dampening effect on the NMR radiation in the active areas of the brain is reduced, so that NMR signals (nuclear magnetic resonance signals) of the water contained in the brain have a higher intensity. Mainly the protons of the hydrogen in the water become excited. DOH consequently has a dampening effect on the NMR radiation. It is reduced as a result of brain activity and consequently, permits a localization of the brain activity if examinations utilizing NMR methods are performed. This action mechanism is known in the field under the name BOLD effect (Blood Oxygen Level Dependent-effect) and results in susceptibility-sensitive magnetic resonance measurements at a field strength of a static magnetic field of for example 1.5 Tesla in up to about 10% changes of the image brightness in the activated brain regions. The endogenic contrast compound DOH may be replaced by other contrast compounds, which cause a change in the susceptibility. With NMR imaging methods, layers or volumes are selected which, with suitable input of high-frequency radiation pulses and the application of magnetic gradient fields, provide a measuring signal, which is digitized and stored in the measuring computer in a two- or three-dimensional field.
Normally, FFT is used, but the use of other transformations might also be used. From the recorded raw data, the desired image information is reconstructed by way of a two- or three-dimensional Fourier transformation.
A reconstructed slice image consists of pixels, a volume data set consists of voxels. A pixel is a two-dimensional image element, for example a square, of which the image is composed. A voxel is a three-dimensional volume element, for example a block which, because of measuring techniques, has no sharp limit. The dimensions of a pixel are in the range of 1 mm2, those of a voxel are in the range of 1 mm3. The geometries and dimensions are variable.
Although the images of the measured slices are only two-dimensional, the term xe2x80x9cvoxelxe2x80x9d is also used for them By a comparison of the measured signals in each pixel with the time-dependent plot model function, a stimulus-specific neural activation can be detected and localized. A stimulus may be for example a somato-sensoric, an acoustic, a visual or an olefactoric stimulus as well as a mental or motor task. The model function or, respectively, model time series describes the expected signal change of the magnetic resonance signal as a result of neural activation. Those may be derived or example by empirical rules from the paradigm of the respective experiments. It is essential to take into consideration a time delay of the model function with respect to the paradigm (slow reaction of the blood flow to neural activation).
It is already known how brain activation can be represented by images, which were derived from nuclear spin tomographic data. The computation and the display of the activation images are even possible on a real-time basis that is an image can be created from a data set before the next set of data is obtained. The time gap is typically 1 to 3 seconds.
The most widely used method of detecting neural activity is the correlation analysis. It is realized by computation of the cross-correlation coefficients between a reference vector, that is a model time series, and the time sequences of the voxels taken into consideration. The real time capability of a correlation algorithm requires that the time needed for the renewed computation of the correlation coefficients is constant for each new set of magnetic resonance data recorded.
A real time correlation analysis is known from Cox, R. W., Jesmanowicz, A., J. S. Magn. Reson, Med., 33 230, 1955, which supports the suppression of low-frequency noise by means of a detrending procedure. In a detrending procedure, one tries to reduce the effect of non-stimulus induced signal changes into the time sequences measured. This means mathematically that the measurement vector, that is, the vector which includes the measured time sequence of a voxel is split into the sum of two orthogonal vectors. The part, which is described by a linear combination of detrending vectors, is discarded. The detrending vectors, which in a mathematical sense form a basis, include the time series by whose weighted sum the complete low frequency noise part of the measured time series is described (xe2x80x9ccompletexe2x80x9d within the limits of the given detrending vectors). By applying detrending to the measuring and reference vectors entering with the correlation computations, the effect of non-stimulus induced signal changes on the correlation images can be reduced.
With the application of detrending, the correlation coefficient xcfx81 is calculated according to the formula (1).                     ρ        =                                                            x                →                            S                        ⁢                                          r                →                            S                                                          "LeftBracketingBar"                                                x                  →                                S                            "RightBracketingBar"                        ⁢                          "LeftBracketingBar"                                                r                  →                                S                            "RightBracketingBar"                                                          (        1        )            
Wherein:
xcfx81=correlation coefficient
{right arrow over (x)}s=measurement vector after detrending
{right arrow over (r)}s=reference vector after detrending
Furthermore, Cox, R. W. presented in Comput. Biomed. Res. 29, 162, 1996, the program AFNI, which is based on the algorithm of the first-mentioned publication. The known algorithm is based on elements of a matrix which is obtained by Cholesky analysis of another matrix which is composed of the scalar products of the measurement and reference vectors as well as the detrending vectors.
The neural activation measured with a nuclear spin tomographic examination follows the stimulation paradigm, that is the protocol, which represents the timing of the stimuli. Usually an on-off period (that is, a stimulus xe2x80x9conxe2x80x9d phase followed by a stimulus xe2x80x9coffxe2x80x9d phase) is repeated several times.
If a correlation analysis takes all measuring time points into consideration, the sensitivity for changes of the signal-response form decreases with each activation cycle. However, the prior art methods do not provide for a uniform sensitivity for the change of the signal response form.
It is therefore the principal object of the present invention to provide a computer and a nuclear magnetic resonance tomograph, which deliver improved analysis results. It is a particular object to provide for a uniform sensitivity for the changes of the signal response form. Also, a computer and a nuclear magnetic resonance tomograph are to be provided which permit the implementation of such a solution. For this purpose, a new detrending procedure for the reduction of undesired low-frequency components in the measured signals is to be developed, whereby the implementation is made possible.
In a computer for evaluating signals from nuclear magnetic resonance tomography and a tomograph using such a computer, the neuronal activity measured during a nuclear magnetic resonance tomographic examination follows a protocol, which describes the timing of the stimuli. During a correlation analysis, the sensitivity of known means to changes in the signal-response form normally decreases with each activation cycle. According to the invention, better evaluation results can be achieved in the correlation analysis by using a constant number of n data values in a correlation analysis employing a sliding-window technique. As the data measurement progresses, the oldest value is continuously disposed of and replaced by the most recent data value in the calculation.
The computer provides for a uniform sensitivity for the changes of the signal response form by utilizing a combination of the equations (1), (4), and (5) below. In a first step, the vectors {right arrow over (xcex3)} and {right arrow over (xcex4)} are calculated in accordance with equations (2) and (3) which vectors are composed of the detrending coefficients xcex3i and respectively, xcex4i ({right arrow over (xcex3)}=(xcex3o, xcex3i . . . )T , {right arrow over (xcex4)}=(xcex4o, xcex4i . . . )T. The equations are obtained by minimizing two square functions. The detrending matrix S includes si as column vectors.
xe2x80x83(STS){right arrow over (xcex3)}=ST{right arrow over (X)}xe2x80x83xe2x80x83(2)
(STS){right arrow over (xcex4)}=ST{right arrow over (r)}xe2x80x83xe2x80x83(3)
wherein:
S=detrending matrix with {right arrow over (s)}i as column vectors
ST=transposed of S (matrix is mirrored along its diagonal)
{right arrow over (x)}=measurement vector; includes the measured time series of a voxel
{right arrow over (r)}=reference vector; includes the model time series
T=symbol for transposed of a matrix,
{right arrow over (xcex3)}=vector which includes the coefficients of the detrending vectors for the calculations of {right arrow over (x)}s.
{right arrow over (xcex4)}=vector which includes the coefficients of the detrending vectors for the calculation of {right arrow over (r)}s.
These linear equation systems are solved mathematically efficiently for example by the method of conjugate gradients, if the reference vector varies spatially, or, otherwise, in the simplest way, by multiplication of the inverse matrix which can be determined by means of the method of conjugate gradients. Other methods may be used which are however less efficient. It is possible to transform the extended correlation equation until only such terms remain which can be iteratively calculated. In this way, the correlation coefficients themselves can be iteratively calculated.
In a second step, the expressions for {right arrow over (x)}s and {right arrow over (r)}s are substituted in accordance with the formula (4) and (5):                                           x            →                    S                =                              x            →                    -                                    ∑                              i                =                0                                            L                -                1                                      ⁢                                          γ                i                            ⁢                                                s                  →                                i                                                                        (        4        )                                                      r            →                    S                =                              r            →                    -                                    ∑                              i                =                0                                            L                -                1                                      ⁢                                          δ                i                            ⁢                                                s                  →                                i                                                                        (        5        )            
wherein:
{right arrow over (x)}=measurement vector; includes the measured time series of a voxel
L=number of detrending vectors
xcex3i=coefficients of the detrending vectors for the calculation of xs 
{right arrow over (r)}=reference vector
xcex4i=coefficients of the detrending vectors for the calculation of rs 
{right arrow over (x)}s=measurement vector after detrending
{right arrow over (r)}s=reference vector after detrending
{right arrow over (S)}i=detrending vectors
S=detrending matrix with Si being column vectors
By combining the formulae (1), (4), and (5), formula (6) is derived which represents the basis for the calculation of the correlation coefficients in accordance with the invention, such. that a uniform sensitivity of the calculated values for the changes of the signal response form is guaranteed. The expression for the correlation obtained therewith can be formulated in terms of sub-expressions, which can be iteratively calculated.
In order to guarantee now a uniformly high sensitivity for changes of the signal response form, the xe2x80x9csliding-Windowxe2x80x9d technique was developed, wherein the correlation is limited in accordance with the invention to the n last magnetic resonance data sets. The number n of data sets included in the calculation can be freely selected. Preferably, the number n is maintained constant during a measuring procedure but, while maintaining the iterative computability of xcfx81, the window width can be increased or reduced by one in each time step, or, with a correspondingly increased calculation effort, the window width can be changed as desired The correlation computation includes the respective newly added data value while the oldest data value contained in the data window of n elements used for the computation is discarded. A data window always consists of n data sets, which may be especially two- or three-dimensional, wherein n may be for example between 10 and 20. The selection of this window width (=number of data employed) is decided upon by the person performing the experiments and depends, for example, on the desired sensitivity for the changes of the measured neural activation and the time between two measurements.
A data value in the sense of the invention is to be understood to be an individual number such as an integer number value or an floating point value measured at a particular time point. One data value is assigned to one voxel. A data set is a collection of data values, which belong to a measuring time. A data set comprises the data values of all voxels or a partial amount of all voxels. A data set may comprise one or several slices.
In the method used by R. W. Cox, the correlation analysis occurs exclusively in a cumulative manner, that is, all data sets recorded from the beginning of the measurements enter the correlation computation. As a result, in this state of the art method, the deviations of the measuring signals from the model functions have a smaller effect on the correlation coefficient the later they occur. In contrast thereto, the sliding-window technique according to the invention provides for a uniform sensitivity for the changes of the signal response form. During the initial fill-up phase, the sliding-window method corresponds to a cumulative correlation. Thereafter, the phase of the xe2x80x9cstationary equilibriumxe2x80x9d follows (that means the number of data values n remains constant, wherein with progressing measurement, the oldest data value is eliminated while the newest data value is added) which, like in the Cox algorithm, requires for each time step, a constant calculation effort.
The algorithm developed for performing the sliding-window procedure is represented by the general formula (6):                     ρ        =                                                            x                →                            ⁢                              r                →                                      -                                          ∑                i                            ⁢                                                δ                  i                                ⁢                                  x                  →                                ⁢                                                      s                    →                                    i                                                                                                        h                1                            ⁢                              h                2                                                                        (        6        )            
wherein:
{right arrow over (x)}=measurement vector
{right arrow over (r)}=reference vector
i=index variable, which extends over N data sets
xcex4i=coefficients of the detrending vectors for the calculation of {right arrow over (r)}s 
{right arrow over (s)}i=detrending vectors
hi and h2 are defined by the formulae (7) and (8):                               h          1                =                                            x              →                        2                    -                      2            ⁢                                          ∑                i                            ⁢                                                γ                  i                                ⁢                                  x                  →                                ⁢                                                      s                    →                                    i                                                              +                      2            ⁢                                          ∑                                  i                   greater than                   j                                            ⁢                                                γ                  i                                ⁢                                  γ                  j                                ⁢                                                      s                    →                                    i                                ⁢                                                      s                    →                                    j                                                              +                                    ∑              i                        ⁢                                          γ                i                2                            ⁢                                                s                  →                                i                2                                                                        (        7        )                                          h          2                =                                            r              →                        2                    -                      2            ⁢                                          ∑                i                            ⁢                                                δ                  i                                ⁢                                  r                  →                                ⁢                                                      s                    →                                    i                                                              +                      2            ⁢                                          ∑                                  i                   greater than                   j                                            ⁢                                                δ                  i                                ⁢                                  δ                  j                                ⁢                                                      s                    →                                    i                                ⁢                                                      s                    →                                    j                                                              +                                    ∑              i                        ⁢                                          δ                i                2                            ⁢                                                s                  →                                i                2                                                                        (        8        )            
wherein:
xcex3i, xcex3j=detrending coefficients for the calculation of xs 
sj=detrending vectors s=detrending matrix which includes the si as column vectors.
In a specialized formulation, xcfx81 is obtained in accordance with formula (9):                     ρ        =                              a            -                          b              ⁢                              c                N                                                                                        d                -                                                      b                    2                                    N                                                      ⁢                                          e                -                                                      c                    2                                    N                                                                                        (        9        )            
wherein:
N=number of measuring points in the measuring window or current number of measuring points (particularly during the initial phase of the correlation computations)
In the equations (10) to (14) a, b, c, d and e are defined as follows:                     a        =                              ∑            i                    ⁢                                    x              i                        ⁢                          r              i                                                          (        10        )                                b        =                                            ∑              i                        ⁢                          r              i                                =                      N            ⁢                          r              _                                                          (        11        )                                c        =                                            ∑              i                        ⁢                          r              i                                =                      N            ⁢                          x              _                                                          (        12        )                                d        =                              ∑            i                    ⁢                      r            i            2                                              (        13        )                                e        =                              ∑            i                    ⁢                      x            i            2                                              (        14        )            
wherein:
{right arrow over (x)}=average value of N measured values
{right arrow over (r)}=average value of N components of the reference vector
ri components of the reference vector
xi components of the measurement vector
The data sets entering the correlation computation for the equations (6) and (9) as well as (1) can be weighted depending on their index number.
The computation time for a correlation image in accordance with the algorithm according to the invention following formula (6) does not depend on the window width, that is the number n (N may indicate any desired number of data values or data sets; if the sliding-window method is used, for the stationary case N=n).
If there is a (time-dependent) shift between measurement and detrending vector, the scalar product between the vectors can no longer be calculated iteratively. As a result, the detrending vectors can no longer be constant but must be moved as sliding-windows over predetermined number series.
For the detrending vectors, the following selection can be made:
It is assumed that the detrending vectors are the vectors with the elements ti, wherein i numbers the detrending vectors and t runs from the smallest to the largest index number of all magnetic resonance data sets taken into consideration for the correlation. These detrending vectors are called polynomial. The series ti are increasing monotonically. This resultsxe2x80x94disregarding small fluctuationsxe2x80x94in a general decrease of the detrending coefficient during the measurement procedure, since they have to compensate for the increase of the series. A direct visualization of the detrending coefficient would therefore be meaningless. Independent of the experiments performed, similar decreases of the detrending coefficient during the measurement would then occur. These coefficients would provide information concerning the time-dependent changes of the low frequency noise components, which would be difficult to interpret.
In the present case, the coefficients can be easily transformed, in accordance with formula (15) into the values xcex3k which one would obtain with a set of constant detrending vectors:                                                         γ              ~                        k                    =                                    ∑                              i                ≥                k                                            L                -                1                                      ⁢                                                            γ                  i                                ⁡                                  (                                                                                    i                                                                                                            k                                                                              )                                            ⁢                              t                0                                  i                  -                  k                                                                    ,                            (        15        )            
wherein:
xcex3k=corrected xcex3-coefficients
L=number of detrending vectors
i=index variable
k=index variable
to=origin of the sliding-window relative to the beginning of the measurement
xcex3i=uncorrected xcex3 coefficient
the index to is the index of the oldest data set in the window (buffer).
The matrix multiplication is utilized for the vector with the detrending coefficients in each voxel for each time step. The computation efforts for the correction of the detrending coefficients are typically smaller than those for the actual correlation. The subsequent linear transformation of the detrending coefficients, which were determined by the systems of linear equation systems, are required for obtaining from the coefficients information in addition to the correlation coefficients. The additional information comprises particulars concerning the strength of the individual low frequency noise components.
For the instance of non-polynomial detrending vectors, appropriate formulae can be derived. By Taylor expansion of the functions on which the non-polynomial detrending vectors are based, the more general case can be reduced to the case of polynomial detrending vectors.
The correlation results depend decidedly on the selection of the reference vector. Starting with a reference vector, which provides a correlation value at least for some image points (pixels), which correlation value is above a certain threshold (for example 0.5), an iterative improvement can be performed. To this end, a vector is generated by component-wise averaging of the time sequences, whose correlation values are above the threshold, and this vector is used as reference vector for a new correlation computation. The mathematical formulation of this procedure is given in formula (16):                               T          i                =                                            Σ                                                ∀                                      pixels                    ⁢                                          (                                              k                        ,                        l                                            )                                                                      ⩓                                                      p                    ⁢                                          (                                              k                        ,                        l                                            )                                                         greater than                   t                                                      ⁢                          x                              i                ,                k                ,                l                                                                        Σ                                                ∀                                      pixels                    ⁢                                          (                                              k                        ,                        l                                            )                                                                      ⩓                                                      p                    ⁢                                          (                                              k                        ,                        l                                            )                                                         greater than                   t                                                      ⁢            1                                              (        16        )            
wherein:
k, l=pixel coordinates
xcfx81=correlation coefficient
t=correlation threshold (may be between xe2x88x921 and 1 and must be exceeded by the correlation values for a consideration of the respective measurement vector in the sum). ri=the ith component of the reference vector xi,k,l=the ith component of the time sequence of the pixel with the coordinates k, l.
With the new reference vector, a new correlation computation is performed. This procedure, which comprises the calculation of the new reference vector and the correlation coefficients and which is considered to be an iteration step, can be repeated several times whereby the reference vector is optimized: Again the time series for the generation of the reference vector are averaged whose respective correlation coefficient is above the predetermined threshold value. The optimizing occurs by adaptation of the reference vector to the time series, which have the highest correlation coefficients. The threshold can be changed from iteration step to iteration step for example in order to limit the selection of the time series to be averaged. The procedure can also be applied to partial areas for which there are individual reference vectors.
Evaluation of the method according to the invention:
In the stationary equilibrium, the computation time with a SUN Ultra Sparc (143 MHz) for a correlation image with a matrix size of 64xc3x9764 pixels and a detrending vector with identical entries is 0.13 s. Together, with the time required for the scaling of the original and the correlation image of 256xc3x97256 pixels by means of bi-linear interpolation a rate of 2 fps (frames per second) was obtained. The implementations of the algorithms according to the invention were compared with the correlation module of the AFNI (Analysis of Functional Neuro Images) software package with respect to speed, which according to the program description is based on the algorithm described in the publication of Cox, R. W., Jesmanowicz A., Hyde, J. S., Magn. Reson. Med. 33, 230, 1995. The software package is described in the publication of Cox R. W. Comput. Biomed. Res. 29, 162, 1996. Comparable computing times were determined. (Because of the missing sliding-window capability, the comparisons have however only limited value.) During the speed test of the real-time software developed by the applicants, the correlation is performed only with one reference vector, which is constant during a measurement. This however is no limitation in principle. The computation was performed also with only one detrending vector with constant elements.
In summary, it can be said that the computer equipped with the algorithm according to the invention delivers with comparable computing times better results.
The way the computer operates in accordance with the invention is not limited to real time measurements. It is also possible to employ measurement values, which are stored in data storage devices.
Also data values and data sets from nuclear magnetic resonance spectroscopy measurements can be evaluated by the computer according to the invention. Consequently, also any nuclear magnetic resonance spectrometer should be included in the term nuclear spin tomograph. In accordance with the invention, the computer includes evaluation means, which operate according to the formula (2), (3) and (6) to (16) or the mathematical equivalent expressions thereof. These means could be chips or respectively, integrated circuits, that is, particularly processors or programmable logic components.
With the computer according to the invention and with the nuclear spin tomograph utilizing this computer, it can also be determined if a patient is susceptible to a stroke so that important decisions can be more competently made before surgeries. In addition, therapy plans and surgery planning for stroke patients can be performed and supported.
Furthermore, the blood supply to all organs can be determined.
The essential steps of the operating method of the computer according to the invention and of the nuclear resonance tomograph are here shortly summarized:
First, the number of detrending vectors L and the window width n are determined. Storage for the required auxiliary variables of the equations (2), (3) and (6) or respectively, (9) is allocated. Also, the correct initialization values are determined. For the equation (6) as well as the equation (9), a distinction between the xe2x80x9cfillupxe2x80x9d mode (cumulative correlation at the beginning of the computation) and the stationary phase has to be made.
In the first case, only numbers are added to the auxiliary variables (in equation (9), a, b, c, d, and e are the auxiliary variables; in equation (6) an auxiliary variable is assigned to each sum), in the last case, additionally values are subtracted.
If the equation (6) is used, the detrending coefficient must be calculated using the formulae (2) and (3). The matrices of the systems of equation (dimensions Lxc3x97L) and the vectors on the right side of the equal sign can be iteratively actualized with each time step. After actualization of the auxiliary variables of the equation (6), the correlation coefficient can be calculated according to the formula (6). The solution of the formula (9) follows correspondingly. In this case, equations (2) and (3) are not used. In the examples, the correlation computation for a single correlation value is demonstrated wherein three polynomial detrending vectors were selected (t{circumflex over ( )}0, t{circumflex over ( )}1, t{circumflex over ( )}2). In the first example, the reference and the measurement vectors of an actual experiment were used; in the second example, they were generated by chance.